Mon. Not. R. Astron. Soc. 000, [TIIT2] (2013) Printed 18 April 2013 (MN JAT^X. style file v2.2) 



The formation of systems with closely spaced low-mass 
planets and the application to Kepler-36 

Sijme-Jan Paardekooper^*, Hanno Rein^, Willy Kley^ 

^ DAMTP, University of Cambridge, Wilberforce Road, Cambridge CBS OWA, United Kingdom 
^Institute for Advanced Study, 1 Einstein Drive, Princeton, NJ 08540 
Cf^ ^ Institut fiir Astronomic & Astrophysik, Univcrsitdt Tiibingcn, Auf der Morgenstelle 10, 72076 Tiibingen, Germany 



CN Draft version 18 April 2013 

< 



6 

c3 



> 

o 
en 

> 

X 



ABSTRACT 

The Kepler-36 system consists of two planets that are spaced unusually close together, 
near the 7:6 mean motion resonance. While it is known that mean motion resonances 
can easily form by convergent migration, Kepler-36 is an extreme case due to the close 
spacing and the relatively high planet masses of 4 and 8 times that of the Earth. 
In this paper, we investigate whether such a system can be obtained by interactions 
with the protoplanetary disc. These discs are thought to be turbulent and exhibit 
density fluctuations which might originate from the magneto-rotational instability. We 
adopt a realistic description for stochastic forces due to these density fluctuations and 
perform both long term hydrodynamical and iV-body simulations. Our results show 
that planets in the Kepler-36 mass range can be naturally assembled into a closely 
spaced planetary system for a wide range of migration parameters in a turbulent disc 
similar to the minimum mass solar nebula. The final orbits of our formation scenarios 
tend to be Lagrange stable, even though large parts of the parameter space are chaotic 
and unstable. 



Key words: planets and satellites: formation 
etary systems: protoplanetary discs 



planetary systems: formation - plan- 



1 INTRODUCTION 

Over the past two decades, more than 800 extrasolar planets 
have been discoverecrl and more than 100 multi-planetary 
systems. The most striking characteristic of this collection of 
planets is the enormous variety, for example in semi-major 



axis, ranging from Hot Jupiters like 51 Peg (Mayor & Queloz 



1995 1, with periods on the order of days, to giant planets on 



orbits of more than 500 years like in HR 8799 | Marois et al, 
2008). Planets have been detected in extreme environments: 



in close binary systems, e.g. 7-Cep ( Neuhauser et al.|2007 l, 
in circumbinary orbits, e.g. Kepler-16 (Doyle et al.||20rT|, 



and around pulsars (Wolszczan & Frail 



19921. A theory of 



planet formation somehow has to be able to explain these 
extreme planetary systems. Or, viewed from a different an- 
gle, extreme planetary systems can be useful testbeds for 
various theories of planet formation. 

A recent addition to the class of extreme planets came 



of dissimilar densities spaced extremely close together. The 
inner planet is a Super Earth, while the outer planet is a 
small Neptune-type planet. Their semi-major axis differ only 
by approximately 1% of an astronomical unit (AU), so that 
the outer planet on closest approach appears twice as large 
as the full moon as viewed from the inner planet. In fact, 
the planets are close to, but just outside, a 7:6 mean motion 
resonance (MMR, Carter et al.||2012 |. For convenience, we 
summarise the main properties of the Kepler-36 system in 
Table El 

The existence of a closely-packed planetary system such 
as Kepler-36 raises some important questions about planet 
formation. How is it possible that systems like Kepler-36 
form in such a compact configuration? While the orbital 
configuration of Kepler-36 is not special in terms of period 
ratio (it is not exactly in the 7:6 MMR), it appears to be spe- 



in the form of Kepler-36 (Carteret al. 2012 1: a pair of planets 
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cial in terms of stability. As found by Deck et al. ( 2012 1, the 
system is very close to unstable regions in parameter space. 
For most of the orbital solutions the Lyapunov timescale is 
of the order of a few hundred years, compared to millions 
of years in our own Solar System ( Hayes|20"07 1 . This means 
the system is chaotic, in the sense that the memory of ini- 
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tial conditions is lost on short timescales. Nevertheless the 



system is stable over much longer timescales (Deck et al 



2012). 



Planet b 



Planet c 



How come that the Kepler-36 planets formed exactly 
on such an island of stability? Their different mean densi- 
ties suggest instead that they have formed in different lo- 
cations, which allowed only the outer planet to capture a 
significant amount of gas from the disc out of which the 
planets formed. Their close proximity to the central star sug- 
gests that both have migrated inward a significant distance. 
Combining both these suggestions, a possible formation sce- 
nario involves convergent migration, a mechanism that is 
known to lead to resonant pairs of planets, sparked by the 
discovery of Gliese 876 ( Marcy et al.|2001a Snellgrove et al 



2001a Lee fc Peale||2002a l. Convergent migration into the 
7:6 MMR could lead to a stable planetary system, and, as 



suggested by Deck et al. (20121, subsequent tidal evolution. 



either by interaction with the star (Papaloizou 2011 Lith- 



wick fc Wu|[20T2l [Batygin fc Morbidelli||2013^ or with the 
remnant disc ( Baruteau fc Papaloizou|2013 1 , could drive the 
planets slightly out of the 7:6 resonance. Furthermore, the 
subsequent evolution could also be driven by evaporation of 
the inner planet's atmosphere ( Owen fc Wu|[2013 l. 

Formation mechanisms involving convergent migration 
have difficulties when the inner planet is massive enough to 
significantly affect the surface density of the disc. In that 
case an outer low-mass planet stops migrating before even 
reaching the 2:1 MMR due to interaction with the density 
waves launched by the inner planet ( |Podlewska-Gaca et al.| 
2012 1. Thus, if both planets are massive, forming resonances 



as widely spaced as the 4:3 MMR poses problems ( Rein et al 
|2012[ ). However, since both planets in the Kepler-36 system 
are in the Super-Earth regime and therefore considered low- 
mass in terms of their interaction with the disc, and we 
can be confident that they do not significantly perturb the 
surface density profile. 

Even so, if the planets formed far apart, how did they 
push through all other resonances on their way to the 7:6? 
This is the question we address in this paper. While it is 
known that for sufficiently high disc masses, low-mass plan- 



ets can be pushed into compact configurations (Papaloizou 
|fc Szuszkiewicz|2005| , the need for a very massive disc brings 
its own problems: the planets would need to have formed 
very early on, possibly during the self-gravitating phase of 
the disc, and they subsequently have to survive the fast mi- 
gration (with a migration timescale of less than a thousand 
years) associated with massive discs. 

In this paper, we explore the effects of a turbulent disc 
on the formation of MMRs for pairs of low-mass planets such 
as the Kepler 36 system. If the stochastic kicks experienced 
by the planets due to turbulent density fluctuations in the 
disc can break some of the early, widely spaced resonances, 
compact configurations may be formed in less massive, more 
realistic discs. These stochastic kicks can be caused by dif- 
ferent physical mechanisms. The most promising one is the 



Magneto-Rotational Instability (MRI, see Sect. 2.1 1 which 



is thought to be active in at least some parts of the proto- 
planetary disc (Nelson & Papaloizou 2004). 



We find that such closely spaced configurations are in- 
deed a very natural outcome of planet disc interactions in 



Mass (lO-^A/,) 
Semi-Major axis (AU) 
Eccentricity 
Period (d) 



-1 nq+u.uyy 
^"^■^-0.081 

U.110Ci_g QQj^g 

< 0.04 
13.83989i«;°««g 



Mean density (g cm-^) 7.46lo;5g 



2 42+0.18 
^•^^-0.14 
n 19QQ+0.0016 

U.lZOO_g QQjg 

< 0.04 
16.23855lH°oi^ 



0.89 



+0.07 
-0.05 



Table 1. Properties of the Kcpler-36 planetary system, where the 
stellar mass is Af* = 1.071 Mq. Data adapted from [Carter et aT] 
l[2012|. 



a standard turbulent protoplanetary disc. The systems that 
form in this scenario closely resemble Kepler-36. Some of 
the orbits are chaotic on short timescales but almost all are 
stable over timescales comparable with the age of the sys- 
tem. These results are remarkably unsensitive to the speed 
of convergent migration and to the strength of stochastic 
forces. 

The paper is organised as follows. We briefiy review the 
physics of planet migration and mean motion resonances 
in Sect. |2] We then describe the results of hydrodynamical 
simulations in Sect. [31 and present the results of A''-body 
integrations in Sect. H] We discuss the stability of the system 
in Sect. [5] We finally summarise, discuss and conclude in 
Sect. H 



2 PLANET MIGRATION AND MEAN 
MOTION RESONANCES 

In this section, we give an overview of the relevant physics 
of planet migration and the formation of mean nrotion res- 
onances. 



2.1 Type I planetary migration 



Ever since [Goldreich fc Tremaine| ( |1980[ ) it has been known 
that embedded planets can exchange angular momentum 
and energy with gaseous discs, leading to changes in their 
orbit. Several types of disc-induced planet migration can be 
distinguished, according to the nature of the interaction with 
the disc ( Ward 1997) . 

Low-mass planets interact with the disc in a linear fash- 
ion, exciting density waves but leaving the overall structure 
of the disc intact. The torque F on the planet resulting from 
these density waves can be calculated semi-analytically for 
an isothermal disc (Tanaka et al.|2002| as 



F = -C|^Epr^f7^. (1) 

Here q denotes the mass of the planet in terms of the mass 
of the central star, h is the aspect ratio of the disc, E is the 
surface density, r is the orbital radius and il is the angu- 
lar velocity. Subscripts p denote that quantities have to be 
evaluated at the location of the planet. The constant C is of 
order unity and depends on the local surface density slope, 
with in general C > leading to a negative torque on the 
planet and therefore inward migration. Allowing for non- 
isothermal effects leads to a different C that now depends 
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very strongly on the local temperature gradient, with the 



possibility of C < and outward migration ( Paardekooper 



in pericenters A-ro = ci7i — 1:172 can be expressed as a linear 
combination of 01 and (^2, and is often used to characterize 
resonant behaviour. 

Several celestial bodies in our own Solar System are in 
resonance. For example the Jovian satellites lo, Europa and 
Ganymede are engaged in a so called 1:2:4 Laplace reso- 
nance, while Neptune and Pluto (and all Plutinos) are lo- 
cated in a 3:2 MMR. However, out of the 8 major planets in 
the Solar System not a single pair is presently in a MMR. 
According to the so called Nice model of the early Solar 



System, this might have been different in the past ( Tsiganis 
et al.|2005| ). 



Extrasolar planetary systems, on the other hand, show 
often evidence for multiple planets in a mean motion res- 
onance. The best studied system in Gliese 876, a system 



Lee fc Pcalc'200i; •2002b; "Snellgrove etal."2001b': 'Nelson fc| 
Papaloizou 2002 : Beauge fc Michtchenko 2003. : Vcras 2007| 



fc Mellema1[2006l [Kley fc Crida|[2008| [Paardekooper et al. 
20101. This migration regime is called Type I. 

For high-mass planets, interaction with the disc be- 
comes non-linear, leading to shocks close to the planet and 
the opening of an annular gap around the orbit of the planet 
I Lin &C Papaloizoul|1986 l. This is known as Type II migra- 
tion. The critical parameter determining wh ether wave-like 
interactions with the disc are linear is q/h'^ (Korycansky & 
Papaloizou|1996 1. The most massive planet of Kepler-36 has 
q/h'^ « 0.2 for the canonical value oi h = 0.05, which means 
that wave-like interactions with the disc are expected to be 
linear. The Kepler-36 planets will therefore be subject to 
Type I migration when embedded in a gaseous disc. 

Since the Type I torque is proportional to the mass 
squared (q^, see Eq. fl|, the migration speed is proportional of massive planets in a 2:1 MMR (e.g. IMarcy et al.||2001b 

to the mass of the planet. Multiple planets embedded in a 
disc, if they are of different mass and migration is directed 
inward, will therefore either migrate away from each other (if 
the inner planet is the most massive) of migrate towards each 
other (if the outer planet is the most massive). The latter 
case of convergent migration is interesting for the formation 
of resonances. 

Type I planet migration is usually studied in laminar 
discs, where a Navier-Stokes viscosity may be included to 
obtain an accretion flow in the disc. In reality, accretion will 
be driven by turbulence in the disc, most likely due to the 
MRI (see |Balbus fc Hawley|1991[ ). The effects of turbulence 
on Type I migration can be profound. Turbulent density 
fluctuations introduce a stochastic component in the torque 
from the disc felt by the planet ( Nelson fc Papaloizou|2004 |. 
In absence of an average migration torque, planets would un- 
dergo a random walk through the disc. If a smooth inward 
migration torque is present, planets would still migrate in- 
ward on average, but if the stochastic component of the 
torque can dominate for long enough, the survival probabil- 



Crida et al. (2008} [Rein fc Papaloizou|2009| ). All these stud- 
ies suggest that a dissipative process that changes the en- 
ergy (semi-major axis) of the orbits is required to bring the 
systems to a resonant configuration. The prime mechanism 
for this type of sculpting of planetary systems is planetary 
migration ( Kley fc Nelson|2012 |. Convergent migration has 
also been shown to be plausible for more closely spaced sys- 



tems such as the 3:2 MMR (Rein et al. 20101, while the 



ity of planets can be increased ( Nelson & Papaloizou 2004 



formation of even more closely spaced 4:3 resonances can 
pose serious problems ( [Rein et al.|2012] ). 

While the protoplanetary disc, and therefore migration 
forces, are present, the planets move within the disc in a 
self-similar fashion (see e.g. Lee fc Peale||2001 1. The plane- 
tary systems formed by such convergent migration process 
typically have both resonant angles 01,2 in small ampli- 
tude libration, a state called apsidal corotation. These plan- 
etary systems are in a stable configuration if the planetary 
masses are not too large. Rarely, resonant systems formed 



Rein fc Papaloizou|[2009l [Adams fc Bloch|2009 l 

It is important to note that there are other possible 
mechanisms that could add a stochastic component to the 
migration forces even in the absence of the MRI. For exam- 
ple, planetesimals and small protoplanets can have a simi- 
lar effect. This makes our discussion slightly more general, 
although the details might vary between different driving 
mechanisms. 

2.2 Mean motion resonances 

When two or more bodies orbit the same central object, 
mean motion resonances (MMR) can occur, where the or- 
bital period ratio of two of the bodies is close to an integer 
ratio. Formally, a system is said to be in a p:q resonance if 
at least one of the resonant angles, 01,2, is librating, i.e. has 
a dynamical range smaller than 27r. The resonant angles are 
defined as 



by migration show signs of an instability ( Pierens fc Nelson 
'2008). Whether the system becomes unstable or not depends 



mainly on the mass ratio of the planets ( Rein et aL][2010 1 
We did not observe such an instability for Kepler-36 in our 
migration scenario. 

Another possibility to lose a resonance after being 
formed by migration is the presence of a stochastic force 
(Rein fc Papaloizou 2009) . In this case the amplitudes of 
both resonance angles 0i,2 undergo a random walk and grow 
proportionally to the square root of time. The diffusion co- 
efficient of this random walk can be estimated analytically 



61,2 = pAi — q\2 + {p — q) ^1,2, 



(2) 



where Ai,2 is the longitude of the inner and outer planet, 
respectively, and ra7i,2 are their pericenters. The difference 



as shown by Rein fc Papaloizou ( 2009 1 or calibrated to nu- 
merical simulations ( [Nelson fc Papaloizou|2004[ [Oishi et al.[ 
[2007|. 

Here, we try to form a system in an extremely close, 
almost orbit crossing, 7:6 resonance. Although the system is 
not formally in resonance, the close proximity suggests that 
it might have been in resonance at some point in its history 
or that the resonance played at least an important role in 
its formation. These results will give us valuable information 
about the possible formation paths of this system and the 
properties of protoplanetary discs during planet formation. 
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3 HYDRODYNAMICAL SIMULATIONS 

We start by describing results from hydrodynamical sim- 
ulations, in which we evolve the gaseous disc dynamically 
together with the planetary system. 



3.1 Numerical method 

We use the publicly available code FARGCJ'^I (Fast Advection 
in Rotating Gaseous Objects, !Massetl2000), which solves the 
locally isothermal equations of hydrodynamics on a cylin- 
drical grid in two dimensions (r, ip) . All quantities (sur- 
face density E, velocity v and pressure p) are treated as 
vertically integrated. As computational domain we choose 
0.2 < r/ro < 1.5 with the full 27r in azimuth. Here, ro is a 
reference radius, which we take to be 1 AU. This domain is 
covered with a grid that has 256 cells in the radial and 768 
cells in the azimuthal direction. 

A surface density profile E — Eo(r/ro)~"'^" was 
adopted, together with a constant kinematic viscosity v — 
lO'^rQ^o- Note that the choice of surface density slope 
implies no viscous evolution of the surface density profile, 
which is convenient when following the evolution of the sys- 
tem for time spans comparable to the viscous time scale at 
the inner edge of the computational domain. In this work, we 
use viscosity primarily for keeping the solution well-behaved 
over long time scales. If we take ro to be 1 AU, the Minimum 
Mass Solar Nebula (MMSN, |Hayashi|[T98T| surface density 
at r = ro equals Eq^mmsn = 2 • 10~ Mt/rg. 

The disc has a constant aspect ratio h = 0.05. As 
mentioned in Sect. |2.1[ this means that the largest planet 



with mass ratio q = Mp/M^ = 2.4 • 10"^ has q/k^ = 0.2, 
which means disc-planet interactions are well inside the lin- 
ear regime as far as wave launched by the planets are con- 
cerned. Corotation torques can be nonlinear for any planet 



mass (Paardekooper & Papaloizou 2009a I, but by taking 



a substantial viscosity together with a locally isothermal 
equation of state we aim at making disc-planet interactions 
as simple as possible for this study. By taking a locally 
isothermal equation of state, we suppress any complications 
from corotation torques due to radial gradients in entropy 



I Paardekooper fc Mellema||2008[ [Baruteau fc MassetpOOS 
Paar dekoope r et al.|2010||Masset fc Casoh|2010[ ). Migration 
is therefore expected to be inward and in the Type I regime. 
Since the two-dimensional equations are obtained by 
averaging vertically, a softening length e of order h must be 
introduced in the gravitational potential of the planets. We 
take e/h = 0.6 for both planets. Decreasing the softening 
length leads to an increase in migration speed and increases 



the relative strength of the corotation torque ( Paardekooper 
fc Papaloizou| |2009a |. A value around e/h — 0.6 gives good 
agreement with three-dimensional simulations ( |Miillcr ct al. 
20T2I). 



3.2 Stochastic forces 

Even with state of the art computing clusters, high com- 
putational costs make it impossible to run fully turbulent 



simulations (for example due to the MRI) for 10* orbits, 
which would be necessary to follow the migration of low- 
mass planets in discs with masses comparable to the MMSN. 
Therefore, we adopt a simplified description of the effects 
of turbulence on the evolution of the orbital parameters in 
terms of stochastic forces ( Rein||2010 1 . Given an autocorre- 
lation time Tc, a discrete first order Markov process can be 
used to generate correlated continuous noise which is then 
added to the force on the planets. The natural scale for such 
a force is the gravitational force per unit mass due to a 
small circular patch of radius rs at a distance \/2rE from 
the planet. Note that this scale, Fo — 7rGS/2, is indepen- 



dent of rs ( Rein|[2010 1 . The correlation time is taken to be 
Tc{r) = Q,{r)~ , and we vary the amplitude F of the corre- 
lated noise from F/Fq = 0.01 to F/Fq = 0.1. These values 
have been calibrated to mimic the forces measured in simu- 
lations of MRI ( Oishi et al.||2007 K However, it is important 
to note that these values are still uncertain. Current MRI 
simulations either simulate the entire disc with low resolu- 
tion, or a local patch with high resolution. The agreement 
between these fundamentally different sets is only marginal. 
The results depend furthermore on many physical proper- 
ties of the disc such as metallicity and ionization fraction. 
As a consequence one has to keep in mind that the value for 
F/Fo might change by more than an order of magnitude. 
The Typ e I migration torque F is proportional to 



{q/hy'Er'^Q^ (Tanaka et al 



2002 



Paardekooper & Pa- 



http : //f argo . in2p3 . f r 



|paloizou]|2009b| l, corresponding to a force per unit mass 
^mig ~ il/h )GS. This means that even the smallest noise 
level we consider {F/Fq — 0.01) gives rise to forces compa- 
rable to disc tides. 



3.3 Results 

3.3.1 Laminar discs 

We first consider the case without stochastic forces. This 
case was also considered in [Papaloizou fc Szuszkiewicz| 
(2005), but without the specific system Kepler-36 in mind. 
We place the outer planet at r/ro = 1-35 and the inner 
planet at r/ro ~ 0.8, which puts them outside of the 2:1 
MMR, and varied the surface density Eq. We have found 
that for all surface densities considered (Eo 5S 2- 10~''M,/rQ, 
the MMSN surface density at 1 AU), the planets migrate 
through the 2:1 resonance, in agreement with the results of 
[Papaloizou fc_Szuszkiewicz| ( |2005| ) and our TV-body simula- 
tions in Sect.B] For Eo ^ 8- 10"'' M« /r^ , the planets end up 
in the 3:2 MMR. Obtaining a closely packed system requires 
a disc that is significantly more massive than the MMSN. 
The results for the highest disc masses, varying from 8 to 32 
times MMSN, are displayed in Fig. [l] 

For Eo = 1.6 • 10"^M,/ro, the planets end up in the 
4:3 MMR (see the bottom panel of Fig. IT]). Disc migration 
is strong enough to push them through the 3:2 MMR around 
f2ot/27r = 850, an event that is recorded as a spike in the 
eccentricity evolution in the middle panel of Fig.fT] After get- 
ting caught in the 4:3 MMR, the eccentricities of the planets 
go up, but only to about 0.01 for the outer, more massive 
planet, and 0.02 for the inner planet. From the top panel of 
Fig. [1] we see that the system continues to migrate inward 
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Figure 1. Evolution of a planetary system with masses similar 
to the Kepler-36 system in a hydrodynamic simulation without 
stochastic forces, for three different values of So that refer to 8, 
16, and 32 times the surface density of the MMSN. Note that 
for T'o = 1 AU, time is in units of years. Top panel: semi-major 
axis; the horizontal dotted line shows the inner edge of the com- 
putational domain. Middle panel: eccentricity, where the largest 
value of e always belongs to the inner, lower mass planet. Bottom 
panel: period ratio; the horizontal dotted lines show first order 
commensurabilities, from 3:2 (top) to 8:7 (bottom). 



at a rate 2 ■ 10~*rono/27r. If we take tq = 1 AU, the asso- 
ciated migration time scale would be only ~ 5000 yrs. This 
is alarmingly short, but we need even higher disc masses, 
and therefore faster migration rate, to come close to the 7:6 
MMR. 

When we further increase the disc mass by a factor of 2 
(Eo = 3.2 • 10"^M*/ro), we find that the planets get locked 
into the 6:5 resonance (bottom panel of Fig. [I|. Again, the 
orbital eccentricities increase, but remain smaller than 0.02, 
in spite of the stronger forcing from the resonance because 
of the stronger eccentricity damping by the disc. Inspection 
of the resonant angles (see Fig. [2| reveals that this configu- 
ration is stable. While (j>i oscillates around a value close to 
zero, (j)2, and therefore <l)2^ 4>i = ^2 — i^i, oscillates around 
a value close to tv, and the oscillations are decreasing in am- 
plitude. We note that the same is true for the 4:3 resonance 
discussed above. 

For So = 6.4- 10~^Af,/ro, the system attains a 8:7 com- 



Figure 2. Evolution of the resonant angle <f>i for the 6:5 resonance 
(top panel), difference between longitudes of periastron (middle 
panel) and period ratio (bottom panel) for a planetary system 
with masses similar to the Kepler-36 system in a hydrodynamic 



simulation without stochastic forces for Eo = 3.2 • 
Note that for ro = 1 AU, time is in units of years. 



10" 



M./r^. 



mensurability (bottom panel of Fig. [TJ . Inspection of the res- 
onant angles reveals that this configuration is stable as well. 
Therefore, it appears that in order to obtain a configuration 
close to the 7:6 MMR, such as Kepler-36, we need a value 
of So that is somewhere in between the highest two disc 
masses we considered. If we take ro = 1 AU, this means that 
we need a disc that is roughly 20 times more massive than 
the MMSN. While such a massive disc is just about gravita- 
tionally stable (Toomre Q ~ 4 for h = 0.05), the associated 
migration time scales are alarmingly short (~ 2000 yrs). 
Moreover, such high disc masses are only expected in the 
very early stages of disc evolution (e.g. Lin fc Pringle|1990 l, 
which means the formation time scales of the planets would 
have to be very short. Furthermore, questions are raised as 
to why none of the planets, while embedded in a very mas- 
sive gas disc, was able to attract a significant amount of gas 
from the disc and thereby become a gas giant planet 

For first-order resonances p-|-l:p, [Papal 
'Szuszkiewicz 



ioizou 



resonances p-|-l:p 
( 2005 1 observed that the dynamics had a 



1 



stochastic character for high values of p, with some reso- 
nances becoming unstable on longer time scales. We have 
not observed such behaviour for our disc and planet masses 
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before the planets reach the inner edge of the grid, which 
may be due to the different mass ratios considered here. In 
any case, we clearly need a very massive disc to come close to 
the 7:6 resonance. However, the exact disc mass needed will 
certainly depend on the adopted surface density and tem- 
perature profile, since these are known to affect the Type 
I migration rates of planets (Paardekooper et al. 20101. It 



should also be kept in mind that unless E oc r~ migration is 
not scale-free. For surface density profiles that are less steep 
than this, migration time scales in terms of the dynamical 
time scale will increase when a planet migrates inward. This 
suggests that closely-packed systems are more easily formed 
at large distances. Note, however, that this effect is rather 
small for a MMSN profile, for which Eo,mmsn oc ^ro/AU. 
Based on our estimate of the necessary surface density Eq 
for a 7:6 resonance to be established, this particular reso- 
nance could have formed in a MMSN disc only at 400 AU. 
For less steep density profiles the situation would be less ex- 
treme. Note that as long as migration remains convergent, 
a given resonance can always be maintained in principle, 
despite migration becoming relatively less strong. 



Papaloizou & Szuszkiewicz ( 2005 1 also noted a depen- 



dence on initial conditions. This is most likely the result of 
incomplete eccentricity damping between resonances. If the 
planets reach & p+ l:p resonance with considerable leftover 
eccentricities obtained in the p:p — 1 resonance, the chances 



of moving through the p+l:p resonance increase ( Kary et al 



1993). This effect becomes more important for higher values 



of p, for which the resonances are spaced closely together. It 
is interesting to note that, as far as this effect is concerned, 
in order to form a closely packed system it can be advan- 
tageous to start the planets farther apart. If we want the 
planets to move through the 6:5 resonance, it is better to 
start the planets outside the 5:4 resonance, so that they ar- 
rive at the 6:5 with some eccentricity, rather than starting 
on circular orbits just outside the 6:5 resonance. Regarding 
this issue, it should be noted that the initial planet orbits in 
FARGO are circular orbits if they could feel only the gravita- 
tional force from the central star. Adding the gravitational 
force due to the disc leads to some initial eccentricity, which 
depends on the disc mass and can be seen at early times in 
the middle panel of Fig. [l] This eccentricity is subsequently 
damped by interaction with the disc, but starting the plan- 
ets too close to a resonance could lead to spurious breaking 
of the resonance due to this initial eccentricity. 



3.3.2 Effect of stochastic torques 

We now let F/Fq ^ 0, and consider first the case of 
So — 1.6 ■ 10~^ Mt/rQ, which was shown above to lead 
to a 4:3 resonance in absence of stochastic forces. In Fig. [3] 
we show the results for four different amplitudes of stochas- 
tic forcing. For F/Fq ^ 0.02, the semi-major axis evolution 
still appears relatively smooth (top panel of Fig. jsl, while 
for higher amplitudes stochastic forces visibly start to affect 
migration. Note that for F/Fo — 0.1, the planets temporar- 
ily swap places around not/2n — 2500. 

Even though the migration rates of the individual plan- 
ets appear unaffected by stochastic forces for F/Fq ^ 0.02, 
they have a profound effect on resonance locking. From the 




1000 2000 3000 4000 5000 

t (2tz/Q.o) 



Figure 3. Evolution of a planetary system witli masses simi- 
lar to the Kepler-36 system in a hydrodynamic simulation with 
So = 1.6 ■ 10~^ with different levels of stochastic forcing. Note 
that for ro = 1 AU, time is in units of years. Top panel: semi- 
major axis; the horizontal dotted line shows the inner edge of 
the computational domain. Middle panel: eccentricity, where the 
largest value of e always belongs to the inner, lower mass planet. 
Bottom panel: period ratio; the horizontal dotted lines show first 
order commensurabilities, from 3:2 (top) to 8:7 (bottom). 



bottom panel of Fig. [31 we see that even the lowest level of 
stochastic forcing moves the planets straight through the 4:3 
resonance. They briefiy settle into 5:4, but this resonance is 
also unstable in the presence of stochastic forcing, and af- 
ter spending some time in 6:5, the planets end up in the 
7:6 resonance before they come too close to the edge of the 
computational domain and we have to stop the simulation. 
Inspection of the resonant angles (see Fig. ffl clearly shows 
the effects of the stochastic kicks experienced by the planets. 
While (/!>i is close to zero on average, and (j>2 close to tt, there 
is no clear libration visible. This is to be expected, since 
stochastic forces continuously try to take the planets out 
of resonance. While the resonance may resist inward kicks, 
if the planets experience a kick that increases their period 
ratio, they are out of resonance, eccentricity damping sets 
in, and the resonant angles start to rotate. Convergent mi- 
gration subsequently takes the planets back into resonance. 
Since the planets are never deep into a particular resonance 
in the presence of stochastic forcing, their eccentricities re- 
main low on average, lower even than when F/Fq = 0. 
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Figure 4. Evolution of the resonant angle 4'1 (for the 7:6 reso- 
nance; top panel), the difference between longitudes of periastron 
(middle panel) and period ratio (bottom panel) for a planetary 
system with masses similar to the Kepler-36 system in a hydrody- 
namic simulation with stochastic forces of amplitude F/Fq = 0.01 
for So = 1.6 • 10~^ M,/r^. Note that for r'o = 1 AU, time is in 
units of years. 



Figure 5. Evolution of the resonant angle {j>i (for the 7:6 res- 
onance; top panel), difference between longitudes of periastron 
(middle panel) and period ratio (bottom panel) for a planetary 
system with masses similar to the Kepler-36 system in a hydrody- 
namic simulation with stochastic forces of amplitude F/Fq = 0.02 
for So = 4.0 ■ 10""* Mt/r^. Note that for ro = 1 AU, time is in 
units of years. 



As soon as the system spends some time in a resonance, 
the longitudes of periastron will be anti-aligned (see the 
middle panel of Fig. p|. If the system gets kicked out of 
resonance, the resonant angles will start to rotate. As soon 
as the next resonance is reached, zu2 — t^i will tend to n 
again. When the resonances are closely spaced, as is the 
case for high values of p, this essentially means that the lon- 
gitudes of periastron remain anti-aligned, which can be seen 
for 9.Qt/{2-K) > 2500 in the middle panel of Fig. [2] 

The presence of stochastic forces allows for closely 
spaced planetary systems in discs comparable to the MMSN. 
This is illustrated in Fig. [S] where we show the resonant an- 
gles in a disc with So = 4.0 • 10"* M«/ro (twice the value 
of the MMSN if ro = 1 AU). Because of the long migration 
time scale in this low-mass disc (f2o''"a/(27r) ~ 10^ — 10^), 
this simulation was done in several steps. As soon as the 
planets came too close to the inner boundary, we stopped 
the simulation, and put the planets back on their original 
position but at the period ratio they had before the restart. 
Care must be taken not to restart too close to a resonance, or 
the system may artificially move through the resonance be- 



cause of the initial conditions (as was discussed at the end 
of Sect. 3.3.11. Time in Fig. [5] is measured since the final 
restart, at which point were just outside the 5:4 resonance. 
After spending ~ 5000 yr (for ro = 1 AU) in the 5:4 reso- 
nance, the system moves into the 6:5, and subsequently the 
7:6. After ~ 21000 yr, the system experiences an outward 
kick, but convergent migration is able to re-establish the 7:6 
resonance after ~ 1000 yr. 

For F/Fo > 0.02 different effects can be seen. Migra- 
tion paths of individual planets become more erratic (see 
top panel of Fig. ^, and consequently migration is no longer 
always convergent. Commensurabilities are no longer readily 
established, but closely spaced systems can still be formed. 
On average, migration is always convergent, pushing the sys- 
tem towards period ratios of order unity, while the strong 
resonances at high p act as a barrier. As a result, for F/Fo — 
0.05, the system hovers between the 4:3 and the 8:7 reso- 
nance for flot/{2Tv) > 2200 in Fig.ls] The strongest stochas- 
tic forcing we consider, F/Fo ~ 0.1, allows the planets to 
switch places, which happens around Q,ot/{2-K) ~ 2500. This 
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hovering between resonances is observed for all surface den- 
sities considered as long as F/Fq > 0.02. 



4 TV-BODY SIMULATIONS 

4.1 Numerical method 

In addition to the hydrodynamical simulations described 
above, we also run A'^-body simulations to study the evo- 
lution of the Kepler-36 system. In comparison to the hydro- 
dynamical simulations discussed in Sect. |3.1[ A'^-body simu- 
lations run much faster and we can explore a wide range of 

parameters. 

We use the freely available rebounet] code ( Rein & 
Liu 20121. The planets are setup in exactly the same way 
as in the hydrodynamical simulations, i.e. just outside the 
2:1 resonance. We add dissipative forces to the equations 
of motion which mimic the interactions of a planet with a 
proto-planetary disc ( Lee fc Peale|[2001 1 . This allows us to 
choose two dimensional parameters for each planet, the mi- 
gration timescale r^ and the eccentricity damping time-scale 
Te, which can be inferred from hydrodynamical simulations 
(Cresswell et al. 20071. In all simulations presented here, 
we use an eccentricity damping timescale which is a fac- 



/\ln6 



tor oi K — 100 shorter than the semi-major axis damping 
timescale. We apply both migration and eccentricity damp- 
ing to the outer planet, but only the eccentricity damping to 
the inner planet. This setup ensures convergent migration. 
Tests have shown that varying K and the precise nature of 
the dissipative forces do not significantly change the results. 
However, simulations without any eccentricity damping on 
the inner planet are more likely to become unstable. This 
is a consequence of angular momentum conservation which 
makes the eccentricity of the undamped inner planet grow 
while the outer planet pushes it inwards. Stochastic forces 
are added in the same way as in the hydrodynamical simula- 



tions (see Sect. 3.2 1. As a normalisation, we use a disc with 
a surface density of Eo = 1.6 ■ 10~^ M«/ro which allows us 
to use the same notation (F/Fo) as in the previous section. 



4.2 Results 

A'^-simulations allow us to confirm and understand the re- 
sults of the hydrodynamic simulations with a much wider 
set of initial conditions. We ran two sets of simulations, one 
with and one without stochastic forces. 

The simulations without stochastic forces but with 
smooth migration are presented in Fig. [6] The plot shows 
a whole set of simulations. The colour represent the pe- 
riod ratio P2/P1 of the planets. Time is evolving upwards. 
Initially planets start outside the 2:1 resonance (as in the 
hydrodynamical simulations). Due to convergent migration, 
the period ratio shrinks. For all migration rates shown here 
{Ta < lO'' yrs) the planets pass through the 2:1 commensu- 
rability. As one can see from this plot, in a smooth migration 
scenario, only extremely fast migration rates of the order of 
1000 yrs are able to form resonances of high p. 




lo-" 10' 

migration timescale [yrs] 

Figure 6. Results from Af-body simulations with smooth mi- 
gration forces. Colour represents period ratio. Time evolves up- 
wards. Stable resonances form for migration rates longer than a 
few thousand years. Which resonance forms is determined by the 
migration rate. 



https : //github . com/hannorein/rebound 



The picture changes dramatically when stochastic 
forces are included. These results are presented in Fig. W\ 
As in Fig. [6J time evolves upwards and color represents pe- 
riod ratio. The four panels correspond to four different mi- 
gration rates, from top left to bottom right, r^ = 2000 yrs, 
10000 yrs, 20000 yrs and 100000 yrs. For low amplitudes 
of the stochastic forces, F/Fo < 0.01, the results are al- 
most identical to those presented in Fig. [6] When stochastic 
forces are included low-p resonances such as 3:2, 4:3 and 5:4 
become unstable and the planets capture into resonances 
of higher p. This can be seen by the dark blue colour on 
the top right of all four panels. Note that even for moderate 
amplitudes of the stochastic forces, F/Fo > 0.01, resonances 
break up and high-p resonances can form. These resonances 
are not necessarily stable. As one can see in the plots, the 
precise resonance in which planets are in can change fre- 
quently (the color in the plots is not constant) when planets 
are in resonances of high p and stochastic migration forces 
are present. However, the system itself remains stable and no 
planets are lost in most cases. Thus, when planets break out 
of a resonance, they re-capture into a (different) resonance 
as the convergent migration forces are still present. 

This opens up a new route of formation for the Kepler- 
36 system. It is important to point out that it is a stochastic 
process. As one can see from the bottom right panel in Fig. IT] 
which of the high-p resonances forms (7:6, 6:5, 5:4, 4:3) is 
highly unpredictable. However, it is important to point out 
that the systems tend to be in a resonance of high p and stay 
stable for a long period. This in itself is a remarkable result 
as the parameter space contains large regions of instability 
and most orbits are chaotic. The interplay of stochastic and 
dissipative forces stabilises the system and prefers a specific 
set of orbital parameters. 
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Figure 7. Results from A^-body simulations with both smooth and stochastic migration forces. Colour represents period ratio. Time 
evolves upwards. Stable resonances form for migration rates longer than a few thousand years. The stochastic forces assume a disc with 
So = 1.6 ■ 10"'' Mt/rQ. For amplitudes F/Fq larger than 0.01, resonances can break up and planets can capture into high-p resonances. 



5 STABILITY 

In this section we briefly discuss the stabihty of the observed 
Kepler-36 system and our orbital solutions. A much more 
detailed analysis of the long term stability and the resonance 
overlap which can lead to chaotic orbits has been performed 



by Deck et al. (20121 



5.1 Measures of instability 

There are different measures that one can use to define 
wether a planetary system is stable or not which may lead 
to different answers. In this paper, we define three different 
timescales that measure how long a system is stable. 

(i) The Lyapunov timescale is a measure of how fast 
two nearby orbital solutions diverge. If the timescale is short, 
the system is called chaotic (i.e. sensitive to initial condi- 
tions). Often a chaotic system tends to be unstable. We will 
later see that this is not the case for the Kepler-36 system. 



See Wisdom (11983) for more details on the numerical algo- 
rithm. 

(ii) The ejection timescale measures the time until at 



least one planet becomes gravitationally unbound from the 
host star and leaves the systems. The Kepler-36 planets can 
only leave the system when they have a close encounter. 
This measure is therefore closely related to the Hill stability 
of the system. However, as Deck et al. (20121 point out, 



Hill stability is sufficient but not necessary for stability. Our 
definition of the ejection timescale is purely based on direct 
measurements in A^-body simulations. 

(iii) Finally, the Lagrange timescale, as for example 



used by Deck et al. (20121, is the time until at least one of 
the planets' semi-major axis changes by more than 10%. The 
planets might get ejected after they drastically changed their 
orbital parameters. However, as we will see, the more likely 
outcome for parameters similar to the Kepler-36 system is 
that the planet stay in a new configuration and remain stable 
from there onwards. 



5.2 Overview of the parameter space 

We plot an overview of the phase space structure of the 
Kepler-36 system in Fig. [8] The x and y axes correspond 
to the outer planet's initial semi-major axis a2 and eccen- 
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Figure 8. Different meassures of stability from long term Af-body simulations. The maximum integration time is 10® yrs. The inner 
planet is initially on a circular orbit with ai = 0.1153 AU. The x and y axes show the initial eccentricity and semi-major axis of the outer 
planet. The left panel shows the Lyapunov timescale of the system. The middle panel shows the time it takes for at least one planet to 
be ejected from the system. The right panel shows the time it takes for at least one planet to change its semi-major axis by more than 
10%. The black line at a2 = 0.1283 AU are solutions fitted to the observational data found by |Carter et al.| | |2012[ |. Green dots represent 
the stable subset identified by [Deck et al.||[2012[l. The purple dot is the final configuration of our hydrodynamical formation model. 



tricity 62, respectively. The inner planet is initially on a cir- 
cular orbit at a semi-major axis of ai = 0.1153 AU. The 
system is coplanar and all angles are chosen from a uni- 
form distribution. We integrate the 10000 initial conditions 
shown in Fig.lslfor 10^ years with the symplectic integrator 
of REBOUND (see Sect. 4.1 1. The color illustrates the differ- 



Planet b Planet c 



ent measures of stability defined in the previous Sect. |5.1| 
The left plot shows the Lyapunov timescale, the middle plot 
shows the time until one planet gets ejected and the right 
plot shows the time until a semi-major axis changes by more 
than 10%. 

One can see that a large fraction of systems is long term 
stable in the sense that no planet gets ejected (blue solutions 
in middle panel). However, the Lyapunov timescales are re- 
markably short, only a few hundred to a few thousand years 
(left panel). This shows that these, although stable, are in- 
deed chaotic ( [Deck et al!]|2012 K 



A smaller subset of the long term stable solutions is 
Lagrange stable, in the sense that their orbital parameters 
do not change significantly (blue solutions in right panel). 
Their Lyapunov timescale is nevertheless very short (left 
panel). Even Lagrange stable orbits are therefore chaotic. 

There are a few minor caveats when comparing the or- 
bital solutions in such a two dimensional plot as the pa- 
rameter space is higher dimensional. In particular, we are 
marginalising over the inner planet's eccentricity, the mass 
ratio, the mutual inclination, the longitudes and the perias- 
trons. We ran additional test which have shown that changes 
in the parameters which are fixed do not qualitatively change 
the structure of the phase space, but might very well shift 
the precise location of stable regions by a small amount. 
Nevertheless, Fig. [8] gives us valuable information about the 
nature of the system. 



Semi-Major axis (AU) 

Eccentricity 

Longitude of periastron (rad) 

True anomaly (rad) 



0.11530 0.12809 

0.014667 0.008180 

1.8586 -1.2668 

-1.2091 1.2825 



Table 2. Orbital parameters derived from a hydrodynamical sim- 
ulation with Eo = 4.0 • 10"'' M,/rl and F/Fo = 0.02 after the 
surface density has been decreased exponentially to zero. 



5.3 Orbital solution from observations 

The best fits to the transit light curves of Kepler-36 are 
shown in Fig. IS] as a black line. These solutions are the out- 
come of a MCMC regression by Carter et al. (20121 which 



does not take into account long term stability. We scaled 
the semi-major axis such that the inner planet's semi-major 
axis is always at ai = 0.1153 AU. This is effectively fix- 
ing the stellar mass and allows us to over-plot these solu- 
tions in our parameter space survey. Note that by doing so 
the semi-major axis of the outer planet is almost constant 
(02 ~ 0.01283 AU) and the set of solutions appears therefore 
as a vertical line. This is a result of Kepler measuring the 
period (and therefore the period ratio) of transits to high 
accuracy. 

We also over-plot the long term stable solutions found 
by Deck et al. (20121 as green dots. These are a subset of 



the Carter et al. (20121 solutions. The stable solutions are 



outside the nominal 7:6 commensurability and very close to 
the boundary where long term stable solutions can exist. 



5.4 Orbital solution from formation model 

We have taken the hydrodynamic simulation defined by 
Eo = 4.0 ■ 10"* M*/rl and F/Fo = 0.02 shortly after the 
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7:6 resonance is reached, and decreased the surface density 
exponentially over a time scale of Q,ot — 2007r (100 years for 
ro = 1 AU). Thus, convergent migration as well as stochastic 
forcing is turned off slowly. The system then settles into an 
orbital configuration close to the 7:6 resonance. We then set 
ro = 0.336467 AU so that the semi-major axis of the inner 
planet is matching the observed value of 0.1153 AU. The re- 
sulting orbital parameters are summarised in Table p] Since 
the effects of the gas disc are turned off, this solution can be 
continued with an A^-body algorithm. Note that the semi- 
major axis ratio is consistent with the confidence interval 
given in 



Carter et al. (2012k (their table S7). 



The purple dot in Fig. |8| corresponds to the FARGO so- 
lution presented in Table 121 This solution is stable for at 
least 10^ years. Remarkably, the Lyapunov timescale of this 
solution is very long (> 10^ yrs). This suggests that the 



solution is not chaotic. Deck et al. (20121 report that only 



a small fraction (< 1%) of their solutions exhibit such long 
Lyapunov timescales. We can confirm this. All simulations in 
our parameter space survey have a short Lyapunov timescale 
(see left panel in Fig.|8|. It is therefore even more remarkable 
that the stochastic migration scenario picks out a solution 
which is not only stable, but also not chaotic. 

We further tested the robustness of this result by a se- 
ries of A*'-body integrations randomly perturbing the orbital 
parameters given in Table[2]by 1%. The results indicate that 
roughly 50% of those solutions are also not chaotic. We are 
therefore confident that our solution has a finite measure 
in parameter space and can act as an attractor in formation 
scenarios such as the one presented here. Note however, that 
although the best fits from the observations and our solution 
coming from a formation scenario are very similar, they are 
not identical. While the semi-major axis ratio is consistent 
with the observations, inspection of the eccentricity param- 



eters given in table S7 from Carter et al. ( 2012 1 reveals that 



the eccentricities of the hydrodynamic solution are slightly 
too low. We tested that, keeping the other orbital param- 
eters fixed, there are still stable, non-chaotic solutions at 
higher eccentricities that are consistent with [Carter et al.| 
( 2012 1. It is possible that for a different value of Eo or F/Fq, 
higher eccentricities can be reached. We have made no at- 
tempt to fine-tune the hydrodynamic outcome to match the 
observations exactly. 



6 CONCLUSIONS 

In this paper, we have studied a possible formation scenario 
for closely spaced low-mass planets, like Kepler-36, through 
convergent migration in a turbulent disc. 

The most important result of our study is that it is in- 
deed possible to form such a closely-packed system by con- 
vergent migration in the presence of stochastic forces. This 
is the first successful formation scenario for such a highly 
compact system of low mass planets. While a scenario in- 
volving smooth migration requires a very large disc mass, 
adding the effect of turbulent density fluctuations through 
stochastic forces onto the planets naturally lead to the for- 
mation of resonances of very high values of p. Hydrodynamic 
and modified A''-body results give very similar results. We 



find that even though the planets get very close together and 
individual resonances may become unstable, the system it- 
self remains surprisingly stable. When slowly removing the 
gaseous disc, the hydrodynamic simulations tend to pick out 
stable, non-chaotic solutions. 

In order to make the problem tractable computation- 
ally, several simplifications had to be made. First of all, we 
adopted a simple prescription for the effects of turbulence 
through stochastic forces acting on the planets. In a real 
disc, the turbulence is likely created by the MRI. Unfortu- 
nately, global MRI simulations spanning 10* orbits or more, 
necessary to study the formation of the Kepler-36 system at 
low disc masses, are not feasible right now. We have tried to 
choose the parameters of our stochastic forcing, amplitude 
and correlation time, to be as close to those of the MRI as 



possible (see Sect. 3.2 I, but it is clearly necessary to work 
towards more realistic turbulent simulations. We leave this 
to future work. 

We have worked in the 2D approximation, using only 
vertically integrated quantities like the surface density. Us- 
ing an appropriate value for the softening parameter in the 
planets' gravitational potential, it is possible to a large ex- 
tent to reproduce the results of 3D hydrodynamical calcula- 
tions as far as migration is concerned ( Miiller et al.|[2012 |. 
Moreover, fully 3D hydrodynamical simulations for more 
than 10* orbits are beyond current computational resources. 
This means that we can only consider a restricted region of 
parameter space, namely coplanar orbits. While the mutual 
inclination in Kepler-36 system can be constrained to be 
smaller than 2.5° ( Carter et al.|2012 l, making the system at 
least almost coplanar, a small relative inclination can affect 
the stability of the system (Deck et al.||2012| ). Such effects 
are probably best studied using the A'^-body approach. 

A formation scenario involving migration needs to 
somehow address the question of how and at what point 
migration is stopped. One possibility for low-mass planets 
such as in the Kepler-36 system is the existence of a 'planet 
trap', where temperature and/or density gradients are such 
that the Type I torque changes sign. It remains to be seen in 
what way trapping the inner planet in such a way changes 
the current picture. 

We have kept the mass of the planets fixed. If the differ- 
ent densities of the planets of Kepler-36 came about through 
a different formation location, it is necessary for the outer 
planet to obtain its gaseous envelope before the planets 
come close together. The early evolution of the system will 
strongly depend on the mass evolution of the outer planet. 
For example, if the solid core of the Neptune-type planet is 
less massive than the inner Super-Earth, convergent migra- 
tion will naturally only commence once the outer planet has 
obtained its envelope. 

We also ignore the effects due to strong X-ray and EUV 
irradiation on the planets. The work of |Owen fc Wu| ( |2013[ ) 
shows that the inner planet in the Kepler-36 system might 
have undergone a phase of evaporation. If the mass-loss was 
significant, it might have changed the orbital configuration 
of the planets after the protoplanetary disc had long disap- 
peared. The mass-loss is stronger for planets closer to the 
star and might therefore provide an alternative explanation 
for the different observed densities in the Kepler-36 system. 
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There are many opportunities for future studies in- 
volving Kepler-36. More advanced simulations simulations 
should allow for accretion onto the planet, model the irra- 
diation from the star and provide a self-consistent stopping 
mechanism for the planets. The orbital solution picked by 
our formation scenario is very close to, but not exactly, the 
observed one, there is the hope to exactly match the pa- 
rameters. This will allow us to make testable predictions for 
poorly constrained orbital parameters, similar to what [Rem] 
eTall poTol) did for HD45364. 
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